function svd_basics()
% SVD_BASICS 奇异值分解基础
% 
% 介绍奇异值分解的基本概念、性质和几何意义

fprintf('=== 奇异值分解基础 ===\n\n');

%% 1. SVD的基本定义
fprintf('1. SVD的基本定义\n');
fprintf('对于任意m×n矩阵A，存在分解 A = U*Σ*V^T\n');
fprintf('其中：U是m×m正交矩阵，V是n×n正交矩阵，Σ是m×n对角矩阵\n\n');

% 示例矩阵
A = [3, 2, 2; 2, 3, -2];
fprintf('矩阵 A (%d×%d) =\n', size(A));
disp(A);

% 计算SVD
[U, S, V] = svd(A);
fprintf('U矩阵 (%d×%d) =\n', size(U));
disp(U);
fprintf('奇异值矩阵 S (%d×%d) =\n', size(S));
disp(S);
fprintf('V矩阵 (%d×%d) =\n', size(V));
disp(V);

% 验证分解
A_reconstructed = U * S * V';
reconstruction_error = norm(A - A_reconstructed, 'fro');
fprintf('重构误差: %e\n', reconstruction_error);

% 奇异值
singular_values = diag(S);
fprintf('奇异值: [%.6f, %.6f]\n', singular_values(1:2));

%% 2. SVD与特征值分解的关系
fprintf('\n2. SVD与特征值分解的关系\n');
fprintf('A^T*A的特征值 = A的奇异值的平方\n');
fprintf('A*A^T的特征值 = A的奇异值的平方\n\n');

% 计算A^T*A和A*A^T
AtA = A' * A;
AAt = A * A';

fprintf('A^T*A =\n'); disp(AtA);
fprintf('A*A^T =\n'); disp(AAt);

% 特征值分解
eigenvals_AtA = eig(AtA);
eigenvals_AAt = eig(AAt);

fprintf('A^T*A的特征值: [%.6f, %.6f, %.6f]\n', sort(eigenvals_AtA, 'descend'));
fprintf('A*A^T的特征值: [%.6f, %.6f]\n', sort(eigenvals_AAt, 'descend'));
fprintf('奇异值的平方: [%.6f, %.6f]\n', singular_values(1:2).^2);

% 验证关系
fprintf('验证 σ² = λ(A^T*A):\n');
for i = 1:min(size(A))
    fprintf('σ_%d² = %.6f, λ_%d(A^T*A) = %.6f, 差异: %e\n', ...
            i, singular_values(i)^2, i, eigenvals_AtA(end-i+1), ...
            abs(singular_values(i)^2 - eigenvals_AtA(end-i+1)));
end

%% 3. 几何解释
fprintf('\n3. 几何解释\n');
fprintf('SVD描述了线性变换的几何结构\n\n');

% 2D示例
A_2d = [3, 1; 1, 2];
[U_2d, S_2d, V_2d] = svd(A_2d);

fprintf('2D变换矩阵 A =\n'); disp(A_2d);
fprintf('奇异值: σ₁ = %.4f, σ₂ = %.4f\n', diag(S_2d));

% 单位圆的变换
theta = linspace(0, 2*pi, 100);
unit_circle = [cos(theta); sin(theta)];
transformed_circle = A_2d * unit_circle;

% 绘制几何解释
figure('Name', 'SVD几何解释');
subplot(2, 2, 1);
plot(unit_circle(1, :), unit_circle(2, :), 'b-', 'LineWidth', 2);
axis equal; grid on; title('原始单位圆');
xlabel('x'); ylabel('y');

subplot(2, 2, 2);
plot(transformed_circle(1, :), transformed_circle(2, :), 'r-', 'LineWidth', 2);
axis equal; grid on; title('变换后的椭圆');
xlabel('x'); ylabel('y');

% 绘制奇异向量
subplot(2, 2, 3);
plot(unit_circle(1, :), unit_circle(2, :), 'b-', 'LineWidth', 1);
hold on;
% 右奇异向量（输入空间的主方向）
for i = 1:2
    v = V_2d(:, i);
    plot([0, v(1)], [0, v(2)], 'g-', 'LineWidth', 3);
    text(v(1)*1.1, v(2)*1.1, sprintf('v_%d', i), 'FontSize', 12);
end
axis equal; grid on; title('右奇异向量 (输入空间)');
xlabel('x'); ylabel('y');

subplot(2, 2, 4);
plot(transformed_circle(1, :), transformed_circle(2, :), 'r-', 'LineWidth', 1);
hold on;
% 左奇异向量（输出空间的主方向）
for i = 1:2
    u = U_2d(:, i) * diag(S_2d)(i);
    plot([0, u(1)], [0, u(2)], 'm-', 'LineWidth', 3);
    text(u(1)*1.1, u(2)*1.1, sprintf('σ_%d u_%d', i, i), 'FontSize', 12);
end
axis equal; grid on; title('左奇异向量 (输出空间)');
xlabel('x'); ylabel('y');

%% 4. 矩阵的秩和零空间
fprintf('\n4. 矩阵的秩和零空间\n');

% 秩亏矩阵示例
A_rank_def = [1, 2, 3; 2, 4, 6; 1, 2, 3.001];
fprintf('秩亏矩阵 A =\n'); disp(A_rank_def);

[U_rank, S_rank, V_rank] = svd(A_rank_def);
s_vals = diag(S_rank);
fprintf('奇异值: [%.6f, %.6f, %.6f]\n', s_vals);

% 数值秩
tol = max(size(A_rank_def)) * eps(max(s_vals));
numerical_rank = sum(s_vals > tol);
fprintf('数值秩: %d (容差: %.2e)\n', numerical_rank, tol);
fprintf('理论秩: %d\n', rank(A_rank_def));

% 零空间和列空间
null_space = V_rank(:, numerical_rank+1:end);
column_space = U_rank(:, 1:numerical_rank);

fprintf('零空间维数: %d\n', size(null_space, 2));
fprintf('列空间维数: %d\n', size(column_space, 2));

% 验证零空间
if ~isempty(null_space)
    null_check = A_rank_def * null_space;
    fprintf('零空间验证 ||A*null_vector||: %e\n', norm(null_check, 'fro'));
end

%% 5. 条件数
fprintf('\n5. 条件数\n');
fprintf('条件数 = 最大奇异值 / 最小非零奇异值\n\n');

matrices_cond = {[1, 0; 0, 1], [1, 0; 0, 0.01], [1, 1; 1, 1.001]};
names_cond = {'单位矩阵', '对角矩阵', '接近奇异'};

for i = 1:length(matrices_cond)
    A_cond = matrices_cond{i};
    [~, S_cond, ~] = svd(A_cond);
    s_cond = diag(S_cond);
    
    fprintf('%s:\n', names_cond{i});
    fprintf('奇异值: [%.6f, %.6f]\n', s_cond);
    
    if s_cond(end) > eps
        cond_svd = s_cond(1) / s_cond(end);
        cond_matlab = cond(A_cond);
        fprintf('SVD条件数: %.2e\n', cond_svd);
        fprintf('MATLAB条件数: %.2e\n', cond_matlab);
        fprintf('差异: %e\n', abs(cond_svd - cond_matlab));
    else
        fprintf('矩阵奇异，条件数无穷大\n');
    end
    fprintf('\n');
end

%% 6. 经济型SVD
fprintf('6. 经济型SVD\n');
fprintf('对于m×n矩阵，当m≠n时可以使用经济型SVD节省存储\n\n');

% 高瘦矩阵
A_tall = randn(10, 4);
fprintf('高瘦矩阵 A (%d×%d)\n', size(A_tall));

% 完整SVD
[U_full, S_full, V_full] = svd(A_tall);
fprintf('完整SVD: U(%d×%d), S(%d×%d), V(%d×%d)\n', ...
        size(U_full), size(S_full), size(V_full));

% 经济型SVD
[U_econ, S_econ, V_econ] = svd(A_tall, 'econ');
fprintf('经济型SVD: U(%d×%d), S(%d×%d), V(%d×%d)\n', ...
        size(U_econ), size(S_econ), size(V_econ));

% 存储节省
storage_full = numel(U_full) + numel(S_full) + numel(V_full);
storage_econ = numel(U_econ) + numel(S_econ) + numel(V_econ);
fprintf('存储节省: %.1f%%\n', (1 - storage_econ/storage_full) * 100);

% 验证重构
A_recon_full = U_full * S_full * V_full';
A_recon_econ = U_econ * S_econ * V_econ';
fprintf('完整SVD重构误差: %e\n', norm(A_tall - A_recon_full, 'fro'));
fprintf('经济型SVD重构误差: %e\n', norm(A_tall - A_recon_econ, 'fro'));

%% 7. SVD的唯一性
fprintf('\n7. SVD的唯一性\n');
fprintf('奇异值是唯一的，但奇异向量在重奇异值时不唯一\n\n');

% 重奇异值矩阵
A_mult = diag([3, 2, 2, 1]);
fprintf('对角矩阵 A =\n'); disp(A_mult);

[U_mult, S_mult, V_mult] = svd(A_mult);
s_mult = diag(S_mult);
fprintf('奇异值: [%.1f, %.1f, %.1f, %.1f]\n', s_mult);

% 识别重奇异值
unique_s = unique(round(s_mult, 10));
for val = unique_s'
    multiplicity = sum(abs(s_mult - val) < 1e-10);
    if multiplicity > 1
        fprintf('奇异值 %.1f 的重数: %d\n', val, multiplicity);
    end
end

%% 8. 伪逆
fprintf('\n8. 伪逆 (Moore-Penrose逆)\n');
fprintf('A^+ = V * Σ^+ * U^T\n\n');

% 非方阵示例
A_pinv = [1, 2; 3, 4; 5, 6];
fprintf('矩阵 A (%d×%d) =\n', size(A_pinv));
disp(A_pinv);

[U_pinv, S_pinv, V_pinv] = svd(A_pinv);
s_pinv = diag(S_pinv);

% 构造伪逆
tol_pinv = max(size(A_pinv)) * eps(max(s_pinv));
S_pinv_inv = zeros(size(S_pinv'));
for i = 1:min(size(S_pinv))
    if s_pinv(i) > tol_pinv
        S_pinv_inv(i, i) = 1 / s_pinv(i);
    end
end

A_pseudoinv = V_pinv * S_pinv_inv * U_pinv';
fprintf('伪逆 A^+ (%d×%d) =\n', size(A_pseudoinv));
disp(A_pseudoinv);

% 与MATLAB内置函数比较
A_pinv_matlab = pinv(A_pinv);
fprintf('与MATLAB pinv的差异: %e\n', norm(A_pseudoinv - A_pinv_matlab, 'fro'));

% 验证伪逆性质
fprintf('伪逆性质验证:\n');
fprintf('A*A^+*A - A: %e\n', norm(A_pinv * A_pseudoinv * A_pinv - A_pinv, 'fro'));
fprintf('A^+*A*A^+ - A^+: %e\n', norm(A_pseudoinv * A_pinv * A_pseudoinv - A_pseudoinv, 'fro'));

%% 9. 最小二乘解
fprintf('\n9. 最小二乘解\n');
fprintf('超定系统 Ax = b 的最小二乘解: x = A^+ * b\n\n');

% 超定系统
A_ls = [1, 1; 1, 2; 1, 3; 1, 4];
b_ls = [1.1; 1.9; 3.2; 3.8];

fprintf('系数矩阵 A (%d×%d) =\n', size(A_ls));
disp(A_ls);
fprintf('右端向量 b = [%.1f; %.1f; %.1f; %.1f]\n', b_ls);

% SVD最小二乘解
A_ls_pinv = pinv(A_ls);
x_svd = A_ls_pinv * b_ls;
fprintf('SVD最小二乘解: x = [%.6f; %.6f]\n', x_svd);

% 与其他方法比较
x_normal = (A_ls' * A_ls) \ (A_ls' * b_ls);  % 正规方程
x_qr = A_ls \ b_ls;  % QR分解

fprintf('正规方程解: x = [%.6f; %.6f]\n', x_normal);
fprintf('QR分解解: x = [%.6f; %.6f]\n', x_qr);

% 残差比较
residual_svd = norm(A_ls * x_svd - b_ls);
residual_normal = norm(A_ls * x_normal - b_ls);
residual_qr = norm(A_ls * x_qr - b_ls);

fprintf('残差比较:\n');
fprintf('SVD: %e\n', residual_svd);
fprintf('正规方程: %e\n', residual_normal);
fprintf('QR分解: %e\n', residual_qr);

%% 10. 数值稳定性
fprintf('\n10. 数值稳定性\n');
fprintf('SVD是数值最稳定的矩阵分解方法之一\n\n');

% 病态矩阵
A_ill = hilb(6);
fprintf('Hilbert矩阵 A (6×6)，条件数: %.2e\n', cond(A_ill));

% 添加噪声
noise_level = 1e-10;
A_noisy = A_ill + noise_level * randn(size(A_ill));
fprintf('噪声水平: %e\n', noise_level);

% 比较不同方法的稳定性
[U_ill, S_ill, V_ill] = svd(A_ill);
[U_noisy, S_noisy, V_noisy] = svd(A_noisy);

s_ill = diag(S_ill);
s_noisy = diag(S_noisy);

fprintf('奇异值变化:\n');
for i = 1:6
    change = abs(s_noisy(i) - s_ill(i)) / s_ill(i);
    fprintf('σ_%d: %.2e → %.2e (相对变化: %.2e)\n', ...
            i, s_ill(i), s_noisy(i), change);
end

% 子空间角度
if min(s_ill) > 1e-12  % 确保非奇异
    subspace_angle = subspace(U_ill, U_noisy);
    fprintf('左奇异向量子空间角度: %.2e 弧度\n', subspace_angle);
end

end